{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第一章 随机模拟和概率\n",
    "\n",
    "这个世界充满了不确定因素，这也是这个世界这么有趣的原因之一。对于不确定事件，我们在理论上建立了概率论作为分析手段，提出了数理统计作为观察、验证和建模手段。而计算机的出现，使得这些理论和模型，不但可以在现实世界中验证，也可以在虚拟的数字空间中实现、验证和发展。我们这门课，主要目的就是向大家简单介绍一下，目前人类在这方面积累的理论和技术。由于时间关系，我们的介绍将是纲要性的，很多细节不可能展开，而是有待大家在后续课程，或者根据自身的兴趣和选择进行自学。\n",
    "\n",
    "## 一个简单的例子\n",
    "\n",
    "**例 1.1** 客车有 $n$ 个座位，乘客排队上车，对号入座。假定第一个乘客的车票不慎丢失，该乘客任选一个座位坐下。后面的乘客持票上车后，如果自己的座位空着，则坐在自己的座位上，否则也任选一个空座位坐下。请计算，最后一个乘客上车后，刚好坐在自己座位上的概率。（这里没有说明乘客总数，我们就当乘客也是 $n$ 人好了。）\n",
    "\n",
    "这道题作为一道概率论的问题，不能算难。不过我们今天看看这样的问题，如果用计算机模拟来解，是否有难度。我们这里采用的代码是 Python，为省事具体环境直接采用 Anaconda Jupter-Notebook(安装参见 https://mirror.tuna.tsinghua.edu.cn/help/anaconda/ )，大家习惯就好。\n",
    "\n",
    "\n",
    "首先导入随机数生成器，我们这里固定随机数种子，使得我们产生的随机序列每一次都是相同的。其中的具体原理我们以后的章节再讨论。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.random.seed(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用一个数组来记录各座位的占据情况，第 $i$ 个元素代表第 $i$ 个座位实际上被手持第几号座位票的乘客占据。初始全部赋 $-1$，表示没有被占据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_car_seats(n):\n",
    "    A = np.ones(n, dtype='int32') * -1\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里我们不妨设第 $i$ 个上车的乘客，手持的是 $i$ 号座位号，他理当坐在 $i$ 号座位上，但是因为之前出现过不靠谱的乘客，因此他的座位有可能被占据，此时他就加入不靠谱系列，也开始乱坐。（这里的“不妨设”是否合理？）我们将第 $i$ 个乘客上车坐下的过程用一个函数来模拟："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boarding(i, car_seats):\n",
    "    n = car_seats.size\n",
    "    if i == 0:   # 第一个乘客\n",
    "        s = np.random.randint(0, n)   # 乱坐\n",
    "        car_seats[s] = i;\n",
    "    elif car_seats[i] == -1:   # 第 i 个乘客，他的位置还空着\n",
    "        car_seats[i] = i       # 坐到自己的位置上   \n",
    "    else:                      # 座位被占了\n",
    "        m = np.random.randint(1, n - i + 1)   # 在剩下的位子中随机找一个   \n",
    "        k = 0\n",
    "        for s in range(0, n):                 # 遍历全部座位\n",
    "            if car_seats[s] == -1:            # 对空位计数\n",
    "                k+=1                          \n",
    "            if k == m:                        # 找到第 m 个空位\n",
    "                car_seats[s] = i              # 占据\n",
    "                break\n",
    "    return car_seats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里关键一步是如何模拟“在剩下的空余座位中随机坐一个”。因为座位是有序的线性排列，而我又不想每次刷新当前空余座位列表，所以我这里事实上在最小 $1$ 到最大剩余座位总数 $n - 1$ 之间先取一个随机数 $s$，然后**从头逐个遍历全部座位**，如果每发现一个空座位，我就计数器加一，直到数到第 $s$ 个空座位（必然存在），就占据这个位置。也就是占据第 $s$ 个**空余**的位置。这样做的效率在最坏情况下，是 $\\Theta(n^2)$ 的。如果 $n$ 的规模巨大，那么我们可以通过适当的数据结构来降低这个复杂度。\n",
    "\n",
    "接下去的模拟流程大家以后会经常看见。我先给出一次实验的模拟，也就是给定 $n$ 个座位和 $n$ 个乘客，开始逐一上车，判定最后一个乘客是否坐在自己的位置上。注意计算机模拟的一个倾向是，在不严重影响效率的情况下，模拟过程尽可能贴近实际事件发生过程。尽管不可避免的，如上面的代码，有时确实需要利用概率等价性做一些转换。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_test(n):\n",
    "    car_seats = init_car_seats(n)\n",
    "    for i in range(0, n):\n",
    "        boarding(i, car_seats)\n",
    "    if car_seats[i] == i:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后是批量实验。计算机模拟不是概率推导，而是实实在在的统计计数。对规模为 $n$ 的问题，我们做 $N$ 次实验，统计其成功次数 $s$，当 $N$ 充分大时，\n",
    "$$\n",
    "s / N \\approx P\\{n\\},\n",
    "$$\n",
    "这个就是大数定律，也是我们计算模拟和一切统计计算的基础。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_test(N, n):\n",
    "    success = 0\n",
    "    for i in range(0, N):\n",
    "        if one_test(n):\n",
    "            success += 1\n",
    "    success_rate = success / N        \n",
    "    print(success_rate)\n",
    "    return success_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们来给一些计算实例："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5034\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.5034"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batch_test(5000, 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从上面的实例，我们严重怀疑，这个事件的发生概率，和 $n$ 无关，而且是常数 $0.5$。这个结论似乎有点违背我们的直观？那么如何来验证这一点？再次强调，以上实验，均不是严格的概率证明，那如何体现其可信度呢？虽然少量的实验结果给我们一些暗示，但那显然是不够的。唯一可行的做法，是做大量的实验，不但让 $N$ 足够大，而且要做充分多组的实验，并且以一定的**统计分析**或**可视化**的方式，将结果的正确性展示出来。也就是我们在物理实验中常用的，控制变量法。比如我们考虑将规模从 $n = 5$ 线性增长到 $n = 500$，然后看一下结果的变化："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.512\n",
      "0.46\n",
      "0.521\n",
      "0.511\n",
      "0.509\n",
      "0.495\n",
      "0.515\n",
      "0.5\n",
      "0.492\n",
      "0.51\n",
      "0.496\n",
      "0.501\n",
      "0.507\n",
      "0.515\n",
      "0.484\n",
      "0.497\n",
      "0.509\n",
      "0.51\n",
      "0.495\n",
      "0.472\n",
      "0.507\n",
      "0.5\n",
      "0.478\n",
      "0.525\n",
      "0.516\n",
      "0.49\n",
      "0.508\n",
      "0.469\n",
      "0.491\n",
      "0.49\n",
      "0.498\n",
      "0.492\n",
      "0.519\n",
      "0.513\n",
      "0.49\n",
      "0.479\n",
      "0.481\n",
      "0.51\n",
      "0.491\n",
      "0.483\n",
      "0.498\n",
      "0.53\n",
      "0.514\n",
      "0.534\n",
      "0.509\n",
      "0.497\n",
      "0.5\n",
      "0.493\n",
      "0.497\n",
      "0.505\n",
      "0.515\n",
      "0.506\n",
      "0.511\n",
      "0.483\n",
      "0.499\n",
      "0.491\n",
      "0.491\n",
      "0.504\n",
      "0.497\n",
      "0.496\n",
      "0.507\n",
      "0.505\n",
      "0.48\n",
      "0.495\n",
      "0.497\n",
      "0.473\n",
      "0.508\n",
      "0.492\n",
      "0.493\n",
      "0.501\n",
      "0.485\n",
      "0.508\n",
      "0.487\n",
      "0.502\n",
      "0.527\n",
      "0.471\n",
      "0.478\n",
      "0.521\n",
      "0.483\n",
      "0.492\n",
      "0.497\n",
      "0.489\n",
      "0.494\n",
      "0.486\n",
      "0.473\n",
      "0.463\n",
      "0.494\n",
      "0.495\n",
      "0.478\n",
      "0.484\n",
      "0.524\n",
      "0.484\n",
      "0.507\n",
      "0.5\n",
      "0.479\n",
      "0.497\n",
      "0.52\n",
      "0.51\n",
      "0.509\n",
      "0.504\n",
      "0.499\n"
     ]
    }
   ],
   "source": [
    "N = 1000   # 重复试验次数\n",
    "T = 100    # 总试验组数\n",
    "C0 = 5     # 试验起始规模\n",
    "cases=np.arange(C0, C0 + T + 1)        # 产生全部试验规模\n",
    "sr=[batch_test(N, i) for i in cases]   # 逐个试验并记录"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里注意 `sr` 的生成手段是 Python 的一个特色，可以省略书写一组循环。光看输出似乎是支持我们的结论的，但是最好有个图像更加直观一些。接下去准备图像输出。这里要用到 `matplot` 类库，以后全部用到的类库都会预先在文件初始导入（这里如果不想看到每一次的 `print` 输出，可以在上面的函数定义中注释掉以提高效率，但是看不见输出会焦虑）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像输出参数调整一看就会："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFmFJREFUeJzt3X+sZOV93/H3d3e5QI2BGJr1ll0MlrdpqBMXcwWsHSnXwUkhiuCPuA2oKYlje2UpKHF+VbitnJRKtvKrplGoyybBP6LGhNpRsrI2IRXhqqjGFmtRXANdZ4vtZc0vm3ixtpZ3vb7f/DEzzLmzZ2bOzJ47c+9z3y/pau+cee453/Oc53zmzDM/NjITSVJZtsy7AElS+wx3SSqQ4S5JBTLcJalAhrskFchwl6QCjQ33iLgnIl6IiC8MuT8i4vci4nBEfD4i3th+mZKkSTS5cv8IcP2I+28Adnd/9gIfOvOyJElnYmy4Z+b/BP5uRJObgI9lx2eACyNiR1sFSpImt62FdVwCPF25fbS77NnBhhGxl87VPeeee+5Vu3btamHz87WyssKWLb50AfZFlX3RZ1/0tdEXX/ziF7+emf9wXLs2wj1qltV+p0Fm7gP2ASwuLubBgwdb2Px8LS8vs7S0NO8y1gX7os++6LMv+troi4j4SpN2bTycHgWql+A7gWdaWK8kaUpthPt+4Nbuu2auBV7KzNOmZCRJszN2WiYiPg4sARdHxFHg14GzADLzvwIHgB8HDgPfAt6+VsVKkpoZG+6ZecuY+xP4+dYqkiSdMV/ClqQCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCtQo3CPi+og4FBGHI+L2mvsvjYgHI+LRiPh8RPx4+6VKkpoaG+4RsRW4C7gBuAK4JSKuGGj274H7MvNK4Gbgv7RdqCSpuSZX7lcDhzPzqcw8CdwL3DTQJoHzu79fADzTXomSpElFZo5uEPE24PrMfGf39r8GrsnM2yptdgB/DXwP8ArgrZn5uZp17QX2Amzfvv2qe++9t639mJvjx49z3nnnzbuMdcG+6LMv+uyLvjb64i1vecvnMnNxXLttDdYVNcsGHxFuAT6Smb8bEXuAP46I12fmyqo/ytwH7ANYXFzMpaWlBptf35aXlylhP9pgX/TZF332Rd8s+6LJtMxRYFfl9k5On3Z5B3AfQGY+DJwDXNxGgZKkyTUJ90eA3RFxeUQs0HnBdP9AmyPAdQAR8f10wv1rbRYqSWpubLhn5ingNuB+4Ek674p5PCLuiIgbu81+BXhXRDwGfBz42Rw3mS9JWjNN5tzJzAPAgYFl76v8/gTw5nZLkyRNy0+oSlKBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSpQo3CPiOsj4lBEHI6I24e0+ZcR8UREPB4Rf9JumZKkSWwb1yAitgJ3AT8KHAUeiYj9mflEpc1u4L3AmzPzGxHxvWtVsCRpvCZX7lcDhzPzqcw8CdwL3DTQ5l3AXZn5DYDMfKHdMiVJkxh75Q5cAjxduX0UuGagzT8GiIj/BWwFfiMz/2pwRRGxF9gLsH37dpaXl6coeX05fvx4EfvRBvuiz77osy/6ZtkXTcI9apZlzXp2A0vATuChiHh9Zh5b9UeZ+4B9AIuLi7m0tDRpvevO8vIyJexHG+yLPvuiz77om2VfNJmWOQrsqtzeCTxT0+YvMvM7mfkl4BCdsJckzUGTcH8E2B0Rl0fEAnAzsH+gzZ8DbwGIiIvpTNM81WahkqTmxoZ7Zp4CbgPuB54E7svMxyPijoi4sdvsfuDFiHgCeBD4tcx8ca2KliSN1mTOncw8ABwYWPa+yu8J/HL3R5I0Z35CVZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnu2nAefhg+8IHOv5LqNfoQkzRrDz8My8uwtAR79qxeft11cPIkLCzAAw+svl/zN+zYabYM9w1iM50wowJ8ebmz/Lvf7fy7vFx+f2wkdceuBBvx/Ft34b4RO3GtzfNqdR7HY1SALy11+qDXF36T7PoyyYPvRjnX1+L8m8W+zy3cn3uus4NtPOVe646a9yAcPGE+9rHZ1DPsKmytt10X4NVj8MADa1fDvI/1Rld37E6cOL3dRppea/vZ4qz2fW7h/tWvdnZw1FPuaoj17j+TOdhpTtxp1z+s3mlUT5itW+HDH4ZTp+rraTOc6o7HRz+69g++e/asDnA4/Ri8973TrXtcnaOO9UYO/jZqb7KOwWO3Z0/n90EbaXqt7WeLdfveW95mdsx1WmbUU+5qiG3dChH1gdZ0kEwb0tOsf1S9o7bV5IQ5cgT+4A/q62n7imBwUMPkJ+S0Ne3Z02/3gQ/Ub7et/e0dgyNHhp90F10E73lPe9sad+K2eaHQRj9Nso7qsRtmI02v1T1gVU36wDm47xddNDo7YLrjP9dwHzyow0JsZaVzf2azOdi6zh4X0r2/GTyJ77yz2SCsrn9UvVWDJ8ydd8KLL9YfxN4J8/DDq6+eq/W0fTVUdwU9bNvDtFHTsCBoekxHnRSDD8rbumfE4EkX0TmuKyvtbGtUQE57oTBMG/Pgo9YxzbOCcYG53gx7wJrkHK6uq7rvo7Jj8NnynXcCXPLqJjVH56vYZ+8Hz35lPvRPr+KC8+vvf+mb8NhjkCtAdP4j10yILfCGN3TaHDsGF154+u+9v+u1veD81esbXMdZZ8Hhw91t0f8PYgO47PLOenvrr9b70jfhuWe/zat3nLNqu4P1vu518J3vnP73XzkCX/5SZXvR2Xi17mF9M6yeun1vYtg6J2l37NgxLuwdhClqGrXuuvuaHtO67fbWd+Lb8OyznWMQwI4dcPY5/WPe5PjU1ZErp/dF9XgPbqs3RgdrquqNx9dcWt8vvX2fdlyMajfsvibrrhsXk5h0XEyyvk59k/191STncDUvhh2fwex49fbVYyEClvI4mQfr/m/rVeZ25b6wMLozLzi/0zlNA7w34L9ypLM8AVbg+ef6fze4vpc7lNUnUvUA9Q563aB67DFYWTmH51/orLuu3lEhc+GFnWVUHq3p1n3s2PD+GaynOlirNYw6EUb157AHo8FtNzmx6o7jV45MHkB1x6Bu3bXHdGAcrGoXq4/39lev3k7v+Izql2PHTh9zmecQW1aH9llnVY53dN5UUL0AePkCo1LT4MlezchxF0C9Gscdg+qDSnU/qmNwcB295YP7PmrcTqLaZ9Xzp3oMqsex6cVM0z5rquk5XJcXo45P9ffnnq9ZfxOZOZefq666Kqf1/vdnbt2aCZ1/3//+/n2f/nTmued2li8sZJ59duf3c8/t3Fe3ji1bMs86q9/u7rs791fbT1LDJO0+/enOsrvv7tc9WOso1f0d9Xej+uXd7x7eF03W12v34IMPTl1r0/4cZdgxHbW/W7d2bg873r3jM+p41PXtli0rL4+l6j73xtZgDT/2Y8NrGlZDdX8jOj/j+m/wGFTrG+ynUefBJON23LgYVeO2bZ1jWTc2JzmOZ9JnTeod1xdnMr4H1w9vXMkGGbvu3ufexKgXY5q++Di4jiZzZbB6bn5hAU6cWGFhYcvQ+edxLxxV5/J+4Acmn4NsOp86al4P+jVG9Ns0XV/TOehR7xLo9eeZvMA27JgOjoPq/i4swK23TvfiYN3bM3vbWlkJTp6ET35y9T6/+GLnnT6Dr5385E/CQw8Nr2nY60PVNyBU5+abvDY0WB/Au94Fl146+gXkaeaZxxn2JoYtW/r7Njg2ofm7yHoG37Qxqs9GvZYweN+4c7i33XF5UWdw/W960/PPNPrDJo8Aa/FzJlfu1UezpldUdVcVTdYxan133535znf+v7F/P+l2Jvm7Nq7cq1eHTZ9BjOqLYTWNumps+oxpmj6rq2faYzJqndXlW7Z8t3YfB+uq1jDN8a722bR/X1ffqCvNSa9Cq1fuTY7PsHFRV2tvfYNX8aNqqtYw7BiMO25Nn2UPbqtJXowDHMwGGbthw72pMz2Jq+oG9aRPOZvWNO0AarrdUX8z6fp6J0JvKmLUiVZddxtTMU21OQ4yx08NVk/itR6DwwzbbpMHlVHjb9KxOe5Bv26fmtY+rKZpLhSGTQdN++A2zdRlE4b7GmjjYDU9MdoIvrYDbZjBWt/97umu/te6zjaNq72Nk3ia7U7abty22ngtYtyDflvjYNJnn4NGvQ43zYNbWxeDg5qG+4acc5+Xpp++G6XpHPmZfshjlh/vHpxPvPXWzpzxJB/O2gjvda6aV+1Nt9vG5wtGvd7Q5INKvRp6rz9A/Zhuqy97NdV96K1Xz6j1N30drmm9dedw3VcxrJkmjwBr8dO7cp/V1eVaWasr917baftmllMeme3NJ5Zira7cm1oPz4oGX39o47WOSbY77nWFYX/bZn2D69s0V+4b6cuD2jLJVUqTq6NhZv3x7j174MSJI+zZ89q13ZAaWQ/Pino13HPPl/m5n3vtyzWsdS2jPgE67lnMmZxzs1jfJOYa7m1/XH6jmMUBXw8nt+ZrnsFSrWEeD/qD+z7LC531Yq7hvpG+PGgjWg8ntzRvm/VCZ67hvlk7XdJsbcYLnbm/W2YzdrokrbUt8y5AktQ+w12SCmS4S1KBDHdJKpDhLkkFMtwlqUCNwj0iro+IQxFxOCJuH9HubRGREbHYXomSpEmNDfeI2ArcBdwAXAHcEhFX1LR7JfALwGfbLlKSNJkmV+5XA4cz86nMPAncC9xU0+4/Ar8FfLvF+iRJU2jyCdVLgKcrt48C11QbRMSVwK7M/FRE/OqwFUXEXmAvwPbt21me9MvQ16Hjx48XsR9tsC/67Is++6Jvln3RJNyjZlm+fGfEFuCDwM+OW1Fm7gP2ASwuLuZSAd8Utry8TAn70Qb7os++6LMv+mbZF02mZY4Cuyq3dwLV/337lcDrgeWI+DJwLbDfF1UlaX6ahPsjwO6IuDwiFoCbgf29OzPzpcy8ODMvy8zLgM8AN2bmwTWpWJI01thwz8xTwG3A/cCTwH2Z+XhE3BERN651gZKkyTX6yt/MPAAcGFj2viFtl868LEnSmfATqpJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK1CjcI+L6iDgUEYcj4vaa+385Ip6IiM9HxAMR8Zr2S5UkNTU23CNiK3AXcANwBXBLRFwx0OxRYDEzfxD4BPBbbRcqSWquyZX71cDhzHwqM08C9wI3VRtk5oOZ+a3uzc8AO9stU5I0iW0N2lwCPF25fRS4ZkT7dwB/WXdHROwF9gJs376d5eXlZlWuY8ePHy9iP9pgX/TZF332Rd8s+6JJuEfNsqxtGPHTwCLww3X3Z+Y+YB/A4uJiLi0tNatyHVteXqaE/WiDfdFnX/TZF32z7Ism4X4U2FW5vRN4ZrBRRLwV+HfAD2fmiXbKkyRNo8mc+yPA7oi4PCIWgJuB/dUGEXElcDdwY2a+0H6ZkqRJjA33zDwF3AbcDzwJ3JeZj0fEHRFxY7fZbwPnAf89Iv53ROwfsjpJ0gw0mZYhMw8ABwaWva/y+1tbrkuSdAb8hKokFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklSgRuEeEddHxKGIOBwRt9fcf3ZE/Gn3/s9GxGVtFypJam5suEfEVuAu4AbgCuCWiLhioNk7gG9k5uuADwK/2XahkqTmmly5Xw0czsynMvMkcC9w00Cbm4CPdn//BHBdRER7ZUqSJrGtQZtLgKcrt48C1wxrk5mnIuIl4CLg69VGEbEX2Nu9eTwiDk1T9DpzMQP7uYnZF332RZ990ddGX7ymSaMm4V53BZ5TtCEz9wH7Gmxzw4iIg5m5OO861gP7os++6LMv+mbZF02mZY4Cuyq3dwLPDGsTEduAC4C/a6NASdLkmoT7I8DuiLg8IhaAm4H9A232Az/T/f1twN9k5mlX7pKk2Rg7LdOdQ78NuB/YCtyTmY9HxB3AwczcD/wR8McRcZjOFfvNa1n0OlPUNNMZsi/67Is++6JvZn0RXmBLUnn8hKokFchwl6QCGe4TiIhdEfFgRDwZEY9HxC92l78qIv5HRPxt99/vmXetsxARWyPi0Yj4VPf25d2vn/jb7tdRLMy7xlmJiAsj4hMR8X+742PPJh4Xv9Q9P74QER+PiHM2y9iIiHsi4oWI+EJlWe04iI7f635ty+cj4o1t1mK4T+YU8CuZ+f3AtcDPd7+K4XbggczcDTzQvb0Z/CLwZOX2bwIf7PbDN+h8LcVm8Z+Bv8rMfwK8gU6/bLpxERGXAL8ALGbm6+m8CeNmNs/Y+Ahw/cCyYePgBmB392cv8KFWK8lMf6b8Af4C+FHgELCju2wHcGjetc1g33d2B+qPAJ+i80G2rwPbuvfvAe6fd50z6ovzgS/RfYNCZflmHBe9T6u/is678T4F/PPNNDaAy4AvjBsHwN3ALXXt2vjxyn1K3W++vBL4LLA9M58F6P77vfOrbGbuBP4NsNK9fRFwLDNPdW8fpXOibwavBb4GfLg7TfWHEfEKNuG4yMyvAr8DHAGeBV4CPsfmHRswfBzUfbVLa/1iuE8hIs4DPgm8JzO/Oe96Zi0ifgJ4ITM/V11c03SzvM92G/BG4EOZeSXw/9kEUzB1uvPJNwGXA/8IeAWd6YdBm2VsjLKm54zhPqGIOItOsP+3zPyz7uLnI2JH9/4dwAvzqm9G3gzcGBFfpvMtoT9C50r+wu7XT0D911SU6ihwNDM/2739CTphv9nGBcBbgS9l5tcy8zvAnwFvYvOODRg+Dpp8tcvUDPcJdL/G+I+AJzPzP1Xuqn79ws/QmYsvVma+NzN3ZuZldF4s+5vM/FfAg3S+fgI2QT/0ZOZzwNMR8X3dRdcBT7DJxkXXEeDaiPgH3fOl1xebcmx0DRsH+4Fbu++auRZ4qTd90wY/oTqBiPgh4CHg/9Cfa/63dObd7wMupTO4/0VmboovTouIJeBXM/MnIuK1dK7kXwU8Cvx0Zp6YZ32zEhH/DPhDYAF4Cng7nYunTTcuIuI/AD9F591ljwLvpDOXXPzYiIiPA0t0vtr3eeDXgT+nZhx0H/x+n867a74FvD0zD7ZWi+EuSeVxWkaSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAL9PdY+6KI57x+CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(C0, C0 + T + 1), sr, 'b.')\n",
    "plt.plot([C0, C0 + T + 1], [0.5, 0.5], 'r-')   # 参考线 0.5\n",
    "axs = plt.gca()\n",
    "plt.axis([C0, C0 + T, 0, 1])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里可以看到，实际结果在 $0.5$ 附近波动，并没有随 $n$ 增长而增长。现在我们有理由猜测，最终结果就是 $0.5$，我们无法仅用计算机证明这一点。但是我们至少可以检验，它是否违背概率论的一般定理，比如大数定律。如果它服从大数定律，那么随着重复试验次数 $N$ 的增长，频率的波动应该收窄："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n",
      "0.5\n",
      "0.5333333333333333\n",
      "0.65\n",
      "0.46\n",
      "0.4666666666666667\n",
      "0.5285714285714286\n",
      "0.5125\n",
      "0.4444444444444444\n",
      "0.63\n",
      "0.4818181818181818\n",
      "0.5\n",
      "0.5230769230769231\n",
      "0.5\n",
      "0.5\n",
      "0.55\n",
      "0.5588235294117647\n",
      "0.4722222222222222\n",
      "0.4473684210526316\n",
      "0.51\n",
      "0.5428571428571428\n",
      "0.4772727272727273\n",
      "0.5608695652173913\n",
      "0.5041666666666667\n",
      "0.496\n",
      "0.5269230769230769\n",
      "0.5555555555555556\n",
      "0.5607142857142857\n",
      "0.5103448275862069\n",
      "0.5233333333333333\n",
      "0.5064516129032258\n",
      "0.56875\n",
      "0.49696969696969695\n",
      "0.4970588235294118\n",
      "0.48\n",
      "0.4861111111111111\n",
      "0.47297297297297297\n",
      "0.5236842105263158\n",
      "0.5102564102564102\n",
      "0.4625\n",
      "0.48048780487804876\n",
      "0.48333333333333334\n",
      "0.4744186046511628\n",
      "0.5136363636363637\n",
      "0.5133333333333333\n",
      "0.558695652173913\n",
      "0.5042553191489362\n",
      "0.5020833333333333\n",
      "0.4530612244897959\n",
      "0.49\n",
      "0.49607843137254903\n",
      "0.5423076923076923\n",
      "0.4528301886792453\n",
      "0.5203703703703704\n",
      "0.5054545454545455\n",
      "0.5160714285714286\n",
      "0.4789473684210526\n",
      "0.4517241379310345\n",
      "0.5050847457627119\n",
      "0.5233333333333333\n",
      "0.5163934426229508\n",
      "0.4629032258064516\n",
      "0.5095238095238095\n",
      "0.478125\n",
      "0.49846153846153846\n",
      "0.509090909090909\n",
      "0.5492537313432836\n",
      "0.5352941176470588\n",
      "0.4855072463768116\n",
      "0.5042857142857143\n",
      "0.5084507042253521\n",
      "0.525\n",
      "0.5506849315068493\n",
      "0.47162162162162163\n",
      "0.504\n",
      "0.5026315789473684\n",
      "0.4857142857142857\n",
      "0.532051282051282\n",
      "0.4936708860759494\n",
      "0.4975\n",
      "0.47901234567901235\n",
      "0.525609756097561\n",
      "0.4783132530120482\n",
      "0.49642857142857144\n",
      "0.48705882352941177\n",
      "0.5290697674418605\n",
      "0.49195402298850577\n",
      "0.5\n",
      "0.48764044943820223\n",
      "0.4777777777777778\n",
      "0.5076923076923077\n",
      "0.5282608695652173\n",
      "0.5010752688172043\n",
      "0.5085106382978724\n",
      "0.4842105263157895\n",
      "0.48125\n",
      "0.4979381443298969\n",
      "0.4846938775510204\n",
      "0.4868686868686869\n",
      "0.496\n",
      "0.5089108910891089\n",
      "0.5215686274509804\n",
      "0.48058252427184467\n",
      "0.5086538461538461\n",
      "0.5028571428571429\n",
      "0.5094339622641509\n",
      "0.49439252336448597\n",
      "0.5111111111111111\n",
      "0.5027522935779817\n",
      "0.49363636363636365\n",
      "0.5054054054054054\n",
      "0.4875\n"
     ]
    }
   ],
   "source": [
    "N = np.arange(10, 2000, 10)\n",
    "n = 10\n",
    "sr=[batch_test(i, n) for i in N]\n",
    "plt.plot(N, sr, 'b.')\n",
    "plt.plot([10, 1990], [0.5, 0.5], 'r-')   # 参考线 0.5\n",
    "axs = plt.gca()\n",
    "plt.axis([10, 1990, 0, 1])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们看到确实有这种现象。更进一步，我们甚至可以直接绘制出 `sr` 的分布情况，这里我们问题规模取 $n$，每组试验重复次数为 $N$，总共做 $T$ 组："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500\n",
    "n = 5\n",
    "T = np.arange(1000)\n",
    "sr=[batch_test(N, n) for i in T]\n",
    "plt.hist(sr, bins = 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们现在可以进行充分多的模拟，从而确信 $P_n = 0.5$ 这个结论是对的。当然这不是一个能让张胖子满意的证明，但是在很多实际问题中，我们已经可以得到我们的结论。同时，实验足够多次后得到的结论，本身就是应用统计学存在的依据。\n",
    "\n",
    "**总结：**\n",
    "1. 计算机模拟并不是概率论的证明，它无法给出确信无疑的结论，但是对大多数实际问题，特别是复杂到概率论推理成为巨大挑战的问题，它像物理学实验那样，提供了指示和验证；\n",
    "2. 概率论和数理统计，是计算机模拟的建模和分析基础。没有扎实的概率与数理统计基础，将很难在实际工作中正确并灵活地运用计算机模拟技术；\n",
    "3. 对计算机模拟结果的解释和分析，具有主观性。我们一方面要充分运用统计学的理论和工具进行仔细分析和验证，至少在一些基本事实上不能和概率论有违背；另一方面，适当的可视化能极大提高数据的实际使用意义。\n",
    "\n",
    "## 概率对直观的违背\n",
    "很多现代数学理论，由于其思想的高度抽象和升华，与我们的日常直观是不一致的。大家可能已经开始淡忘极限理论带来的不适感，那是因为进过长期的训练，我们可以改进我们思维模式，从而提升自己的思维的强度和深度。这也是数学改进智商的本质原因之一。在众多数学分支中，概率论对直观的违背是最严重的一个。而通过计算机模拟，有助于我们认识我们自己思维的局限，减轻概率论学习在理解上的困难。或者，这也预示着有这种可能性，对于没有掌握高深概率论知识的人，有可能在具备的基础概率理论水平之后，通过计算机模拟提升自己对随机问题的实际分析能力。这也是我在之前暗示，当我们做计算机模拟时，尽可能少代入自己对模型的概率论理解（因为那可能是错的！），而尽可能照搬现实世界中事件的实际发生过程。\n",
    "\n",
    "这里我们列举一些可能造成概率论和直观违背的例子：\n",
    "1. 书上1.2.1提到的 Monty Hall Problem. 这个例子在电影决胜21点中也出现过；\n",
    "2. 如果没有计划生育政策，且每一个家庭都执行必须有男孩的策略，也就是不断地生，直到出现男孩为止。可以想象，或正如你可能在生活中遇到的那样，会出现7，8个姐妹加一个弟弟的家庭组合。请问在这种情况下，长期以往是否会出现男女性别比例失衡？如果在这样一种落后的社会观念体系中，允许对胎儿进行性别检测，并允许无理由堕胎，请问是否会引起男女性别比例失衡？\n",
    "3. 为什么你不可能中彩票？\n",
    "\n",
    "以上问题，既可以用概率论证明，也可以用计算机模拟验证。对于没有受过足够训练的人，会倾向于更接受后者的结果。万幸，在没有错误的情况下，二者的结论总是一致的。\n",
    "\n",
    "## 动态模拟\n",
    "在很多实际场合，我们关心的一个事件是不停动态变化的。比如各种金融指标，比如战争态势，等等。这些复杂的动态问题的理论模型，往往要涉及概率微分方程等艰深的理论，而近似计算模型的偏差往往也很大。但计算机模拟，在投入足够多计算资源的前提下，却有可能给出非常接近实际情况的结果。至少能提高算力充足方在综合博弈中的胜率。我们会在一些专题分析中纲要性地介绍它们的理论的应用前沿。\n",
    "\n",
    "## 准备工作\n",
    "很明显，我们需要以下的学习准备：\n",
    "1. 概率论基础，对于没有学习过概率论的同学，建议仔细阅读本书第二章；\n",
    "2. 计算机编程基础。我们采用 Python 语言来实现我们的全部模拟。Python语言本身非常接近伪代码和自然语言。我们会在整个授课过程中不断掌握。学习一门语言的唯一途径是去使用它；"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
